/*
Copyright (c) 2009-2010, Warwick Warp Limited
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of the Warwick Warp Limited nor the
      names of its contributors may be used to endorse or promote products
      derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL WARWICK WARP LIMITED BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#ifndef WARRAY_VMATH_H_INCLUDED
#define WARRAY_VMATH_H_INCLUDED

#include <functional>
#include <algorithm>
#include <cmath>
#include <complex>
#include "warray_sys.h"
#include "warray_help.h"

namespace warp {
namespace vm {

template<class T>
void add(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,std::plus<T>());
}

template<class T>
void add_const(size_t n, const T *a, const T& b, T *c)
{
	std::transform(a,a+n,c,std::bind2nd(std::plus<T>(),b));
}

template<class T>
void sub(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,std::minus<T>());
}

template<class T>
void sub_const(size_t n, const T *a, const T& b, T *c)
{
	std::transform(a,a+n,c,std::bind2nd(std::minus<T>(),b));
}

template<class T>
void neg(size_t n, const T *a, T *c)
{
	std::transform(a,a+n,c,std::negate<T>());
}

template<class T>
void mul(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,std::multiplies<T>());
}

template<class T>
void mul_const(size_t n, const T *a, const T& b, T *c)
{
	std::transform(a,a+n,c,std::bind2nd(std::multiplies<T>(),b));
}

template<class T>
void mulconj(size_t n, const std::complex<T> *a, const std::complex<T> *b, std::complex<T> *c)
{
	std::transform(a,a+n,b,c,warp::detail::multiplies_conj<T>());
}

template<class T>
void div(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,std::divides<T>());
}

template<class T>
void div_const(size_t n, const T *a, const T& b, T *c)
{
	std::transform(a,a+n,c,std::bind2nd(std::divides<T>(),b));
}

template<class T>
void square(size_t n, const T *a, T *b)
{
	mul(n,a,a,b);
}

template<class T>
void conj(size_t n, const std::complex<T> *a, std::complex<T> *b)
{
	std::transform(a,a+n,b,warp::detail::choose_conj<T>());
}

template<class T>
void abs(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_abs<T>());
}

template<class T>
void abs(size_t n, const std::complex<T> *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_abs<std::complex<T> >());
}

template<class T>
void inv(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::inv<T>());
}

template<class T>
void sqrt(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_sqrt<T>());
}

template<class T>
void invsqrt(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::invsqrt<T>());
}

template<class T>
void pow(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,warp::detail::choose_pow<T>());
}

template<class T>
void pow(size_t n, const T *a, const T& b, T *c)
{
	std::transform(a,a+n,c,std::bind2nd(warp::detail::choose_pow<T>(),b));
}

template<class T>
void hypot(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,warp::detail::hypot<T>());
}

template<class T>
void exp(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_exp<T>());
}

template<class T>
void ln(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_log<T>());
}

template<class T>
void log10(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_log10<T>());
}

template<class T>
void sin(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_sin<T>());
}

template<class T>
void cos(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_cos<T>());
}

template<class T>
void sincos(size_t n, const T *a, T *b, T *c)
{
	detail::alias_check(b,1,n,c,n,1);
	if(a!=b)
	{
		sin(n,a,b);
		cos(n,a,c);
	}
	else
	{
		cos(n,a,c);
		sin(n,a,b);
	}
}

template<class T>
void cis(size_t n, const T *a, std::complex<T> *b)
{
	std::transform(a,a+n,b,warp::detail::cis<T>());
}

template<class T>
void tan(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_tan<T>());
}

template<class T>
void acos(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_acos<T>());
}

template<class T>
void asin(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_asin<T>());
}

template<class T>
void atan(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,warp::detail::choose_atan<T>());
}

template<class T>
void atan2(size_t n, const T *a, const T *b, T *c)
{
	std::transform(a,a+n,b,c,std::ptr_fun<T,T,T>(std::atan2));
}

template<class T>
void sinh(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,detail::choose_sinh<T>());
}

template<class T>
void cosh(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,detail::choose_cosh<T>());
}

template<class T>
void tanh(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,detail::choose_tanh<T>());
}

template<class T>
void asinh(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,detail::choose_asinh<T>());
}

template<class T>
void acosh(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,detail::choose_acosh<T>());
}

template<class T>
void atanh(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,detail::choose_atanh<T>());
}

template<class T>
void floor(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,std::ptr_fun<T,T>(std::floor));
}

template<class T>
void ceil(size_t n, const T *a, T *b)
{
	std::transform(a,a+n,b,std::ptr_fun<T,T>(std::ceil));
}

// Specialisations
// Please implement both MKL and Accelerate versions if you add to these and ensure 
// that the specialised types are covered by a unit test.
#if defined(USE_MKL) || defined(USE_ACCELERATE)
template<> void add<float>(size_t n, const float *a, const float *b, float *c);
template<> void sub<float>(size_t n, const float *a, const float *b, float *c);
template<> void mul<float>(size_t n, const float *a, const float *b, float *c);
template<> void div<float>(size_t n, const float *a, const float *b, float *c);
template<> void add<double>(size_t n, const double *a, const double *b, double *c);
template<> void sub<double>(size_t n, const double *a, const double *b, double *c);
template<> void mul<double>(size_t n, const double *a, const double *b, double *c);
template<> void div<double>(size_t n, const double *a, const double *b, double *c);
template<> void add< std::complex<float> >(size_t n, const  std::complex<float>  *a, const  std::complex<float>  *b,  std::complex<float>  *c);
template<> void sub< std::complex<float> >(size_t n, const  std::complex<float>  *a, const  std::complex<float>  *b,  std::complex<float>  *c);
template<> void mul< std::complex<float> >(size_t n, const  std::complex<float>  *a, const  std::complex<float>  *b,  std::complex<float>  *c);
template<> void mulconj<float>(size_t n, const  std::complex<float>  *a, const  std::complex<float>  *b,  std::complex<float>  *c);
template<> void div< std::complex<float> >(size_t n, const  std::complex<float>  *a, const  std::complex<float>  *b,  std::complex<float>  *c);
template<> void add< std::complex<double> >(size_t n, const  std::complex<double>  *a, const  std::complex<double>  *b,  std::complex<double>  *c);
template<> void sub< std::complex<double> >(size_t n, const  std::complex<double>  *a, const  std::complex<double>  *b,  std::complex<double>  *c);
template<> void mul< std::complex<double> >(size_t n, const  std::complex<double>  *a, const  std::complex<double>  *b,  std::complex<double>  *c);
template<> void mulconj<double>(size_t n, const  std::complex<double>  *a, const  std::complex<double>  *b,  std::complex<double>  *c);
template<> void div< std::complex<double> >(size_t n, const  std::complex<double>  *a, const  std::complex<double>  *b,  std::complex<double>  *c);
template<> void square<float>(size_t n, const float *a, float *b);
template<> void square<double>(size_t n, const double *a, double *b);
template<> void conj<float>(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void conj<double>(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void abs<float>(size_t n, const float *a, float *b);
template<> void abs<double>(size_t n, const double *a, double *b);
template<> void abs<float>(size_t n, const std::complex<float> *a, float *b);
template<> void abs<double>(size_t n, const std::complex<double> *a, double *b);
template<> void inv<float>(size_t n, const float *a, float *b);
template<> void inv<double>(size_t n, const double *a, double *b);
template<> void sqrt<float>(size_t n, const float *a, float *b);
template<> void sqrt<double>(size_t n, const double *a, double *b);
template<> void sqrt< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void sqrt< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void invsqrt<float>(size_t n, const float *a, float *b);
template<> void invsqrt<double>(size_t n, const double *a, double *b);
template<> void pow<float>(size_t n, const float *a, const float *b, float *c);
template<> void pow<double>(size_t n, const double *a, const double *b, double *c);
template<> void pow< std::complex<float> >(size_t n, const std::complex<float> *a, const std::complex<float> *b, std::complex<float> *c);
template<> void pow< std::complex<double> >(size_t n, const std::complex<double> *a, const std::complex<double> *b, std::complex<double> *c);
template<> void pow<float>(size_t n, const float *a, const float& b, float *c);
template<> void pow<double>(size_t n, const double *a, const double& b, double *c);
template<> void pow< std::complex<float> >(size_t n, const std::complex<float> *a, const std::complex<float>& b, std::complex<float> *c);
template<> void pow< std::complex<double> >(size_t n, const std::complex<double> *a, const std::complex<double>& b, std::complex<double> *c);
template<> void hypot<float>(size_t n, const float *a, const float *b, float *c);
template<> void hypot<double>(size_t n, const double *a, const double *b, double *c);
template<> void exp<float>(size_t n, const float *a, float *b);
template<> void exp<double>(size_t n, const double *a, double *b);
template<> void exp< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void exp< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void ln<float>(size_t n, const float *a, float *b);
template<> void ln<double>(size_t n, const double *a, double *b);
template<> void ln< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void ln< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void log10<float>(size_t n, const float *a, float *b);
template<> void log10<double>(size_t n, const double *a, double *b);
template<> void log10< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void log10< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void sin<float>(size_t n, const float *a, float *b);
template<> void sin<double>(size_t n, const double *a, double *b);
template<> void sin< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void sin< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void cos<float>(size_t n, const float *a, float *b);
template<> void cos<double>(size_t n, const double *a, double *b);
template<> void cos< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void cos< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void sincos<float>(size_t n, const float *a, float *b, float *c);
template<> void sincos<double>(size_t n, const double *a, double *b, double *c);
template<> void cis<float>(size_t n, const float *a, std::complex<float> *b);
template<> void cis<double>(size_t n, const double *a, std::complex<double> *b);
template<> void tan<float>(size_t n, const float *a, float *b);
template<> void tan<double>(size_t n, const double *a, double *b);
template<> void tan< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void tan< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void acos<float>(size_t n, const float *a, float *b);
template<> void acos<double>(size_t n, const double *a, double *b);
template<> void acos< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void acos< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void asin<float>(size_t n, const float *a, float *b);
template<> void asin<double>(size_t n, const double *a, double *b);
template<> void asin< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void asin< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void atan<float>(size_t n, const float *a, float *b);
template<> void atan<double>(size_t n, const double *a, double *b);
template<> void atan< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void atan< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void atan2<float>(size_t n, const float *a, const float *b, float *c);
template<> void atan2<double>(size_t n, const double *a, const double *b, double *c);
template<> void sinh<float>(size_t n, const float *a, float *b);
template<> void sinh<double>(size_t n, const double *a, double *b);
template<> void sinh< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void sinh< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void cosh<float>(size_t n, const float *a, float *b);
template<> void cosh<double>(size_t n, const double *a, double *b);
template<> void cosh< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void cosh< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void tanh<float>(size_t n, const float *a, float *b);
template<> void tanh<double>(size_t n, const double *a, double *b);
template<> void tanh< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void tanh< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void asinh<float>(size_t n, const float *a, float *b);
template<> void asinh<double>(size_t n, const double *a, double *b);
template<> void asinh< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void asinh< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void acosh<float>(size_t n, const float *a, float *b);
template<> void acosh<double>(size_t n, const double *a, double *b);
template<> void acosh< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void acosh< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void atanh<float>(size_t n, const float *a, float *b);
template<> void atanh<double>(size_t n, const double *a, double *b);
template<> void atanh< std::complex<float> >(size_t n, const std::complex<float> *a, std::complex<float> *b);
template<> void atanh< std::complex<double> >(size_t n, const std::complex<double> *a, std::complex<double> *b);
template<> void floor<float>(size_t n, const float *a, float *b);
template<> void floor<double>(size_t n, const double *a, double *b);
template<> void ceil<float>(size_t n, const float *a, float *b);
template<> void ceil<double>(size_t n, const double *a, double *b);
#endif

} // namespace vm
} // namespace warp

#endif
